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Abstract 

The equilibrium phase diagram of a 1:1 symmetrical mixture composed of oppositely 
charged colloids is calculated using Monte Carlo simulations. We model the system 
by the DLVO effective interaction potential. The phase diagram is similar to that of 
its atomic analog (the ionic fluid), where a liquid-gas first order transition emerges 
in the low T — p regions being stable with respect to crystallization. As in the 
ionic fluids, we have found two different crystals: at high T the fluid crystallizes 
in a FCC lattice, whereas at low T, the liquid coexists with a BCC crystal. The 
region of gas-liquid stability is observed to be narrower as the interaction range is 
diminished. 



1 Introduction 

The application of statistical mechanics concepts from atomic systems has 
boosted the understanding of more complex systems in the last decade, in 
particular of phase diagram of macromolecules in solution; i.e. the equilibrium 
phases arise from the competition of the energetic terms and the entropic ones 
[1]. Taking colloids as most simple macromolecules, without internal degrees 
of freedom, they still present some advantages: first, the interaction between 
the particles are easily tailored by external parameters (such as addition of 
salt or polymers), and secondly, the much larger length and time scales in col- 
loids make experiments easier to handle. Additionally, in mixtures of colloids 
different interactions between provoke correlations that may introduce new 
phases. 

For monocomponent colloidal systems, the phase diagrams are quite well un- 
derstood, although not yet completely [1]. An attraction term in the inter- 
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action potential induces fluid-fluid coexistence, whose critical point moves to 
lower temperature as the range of the interaction is decreased. At a certain 
(short) range, the fluid-fluid coexistence becomes metastable with respect to 
crystallization [1,2,3,4]. The crystal phases are also affected, appearing an 
isostructural FCC-FCC transition as the range of the attraction narrows. 

Recently, we studied the colloidal analog of the ionic fluid [5] (Restricted Prim- 
itive Model, RPM), i.e. a symmetrical mixture of oppositely charged colloids. 
As in the RPM, the colloidal mixture presents a gas-liquid separation in the 
low T - low p regions where high correlations between oppositely charged 
particles are found [6]. Furthermore, the critical temperature evolves non- 
monotically as the range of the interaction (salt concentration in the medium) 
decreases [7] . Experimental works with this system have focused on the crystal 
phases of such colloidal mixtures [8,9]. Due to the high correlation between 
opposite colloids, different superlattice crystals are found, with cesium chlo- 
ride structure for symmetrical systems (charge and size) at high attraction 
strengths and disordered-FCC for low ones, when the charge correlations are 
weak. Indeed, these structures are observed in the RPM [10], and we seek 
them in the simulations in this work. 

In this paper, our aim is to study the stability of the liquid for several interac- 
tion ranges, i.e. whether the freezing line crosses the liquid branch. We need 
for this to determine accurately the gas-liquid coexistence line, focusing on the 
finite size effects (on passing, we also show that this transition belongs to the 
3D Ising criticallity class). On the other hand, the freezing line is estimated 
by melting a crystallyte, whose structure is also studied. The results indicate 
the existence of a stable liquid in a narrow region of temperatures for all of 
the ranges studied, and that the crystals have CsCl structure, in agreement 
with experimental results and with the RPM. 

The paper is organized as follows: in Section II, we present the model and the 
simulations details are explained. Section 111 is divided in three subsections, 
where we study the liquid-gas transition in detail, the structure of the crys- 
tals, finally, the freezing line and the stability of the liquid-gas transition. We 
conclude in Section IV presenting the main conclusions. 

2 Model and Simulations 

2.1 Model 

Our system is composed of a 1:1 binary mixture of spherical colloidal par- 
ticles of equal diameter; A^/2 bearing a surface potential +0 and A^/2 with 
— 0. The mixture is immersed in a continuous medium characterized by its 
dielectric constant, £, in presence of an electrolyte. The electrostatic interac- 
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tions are modeled using the well-known DLVO effective potential between the 
colloidal particles [11] and the dispersion forces are not considered. Thus the 
total potential is given by: 



oo, 



r < a 



U{r) = < 



(1) 



7r£(70i02 • exp {—K{r — cr)} , 



r > a 



where a is the hard core diameter of the particles, 0i and 02 are the surface 
potentials, and k, is the inverse Debye length, which depends on the ionic 
concentration. If the charges of the particles are similar, the interaction in (1) 
is repulsive, whereas an attraction results if the particles carry opposite surface 
charges. This system is, thus, the colloidal analogue of the widely studied 
Restricted Primitive Model (RPM). This is obviously an approximation to the 
whole problem, where the ions are not simulated. However, the global problem, 
where the small ions are considered, requires a big amount of computational 
time and the coexistece lines cannot be easily obtained [12]. The results here 
are, thus, a first approximation to the real problem. 

Hereafter, reduced units will be used: a = 1, U* = U / {Tieacj)^) and so T* = 
ksT/ {7iea(j)^) defines the reduced temperature, n* = n — /c^TlogA'^ (A is 
the thermal de Broghe wavelength) is the reduced chemical potential and the 
density is defined as p* — Na^/V (where N is the number of particles, and V 
is the volume of the system) . 

2.2 Finite Size Scaling of Bruce and Wilding 

In the neighborhood of the critical point, the coexistence curve obtained from 
simulations depends strongly on the system size due to the long range of the 
correlations. This fact can be used to locate the critical point using the ap- 
proach developed by Bruce and Wilding (BW) [13]. This method is based on 
the asymmetry of the density distribution, reflecting the absence of particle- 
hole symmetry in off-lattice models. Thus, the order parameter of the tran- 
sition is not the density, but a mixture between the density and the energy: 
M ~ p + su (where s is a system-dependent field mixing parameter). Precisely 
at criticality, the distribution of M in a large enough system, with box size 
L, takes on the universal form: Pl{M) = P*{x = ai^M— < M >)), where 
P* is an universal distribution function for each universahty class, < M > is 
evaluated at critical conditions and ~ L^/'^. Thus, using GCMC simula- 
tions, critical parameters for different box sizes are calculated, and applying 
the corresponding scaling laws, the critical parameters for an infinite system 
can be obtained (Tc(L) — Tc{oo) ~ L-(^-^)/'^^ where u is the critical exponent 
associated to the correlation length and 9 is the universal correction to scahng 
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exponent). This method can be improved inchiding the pressure as a scahng 
field, but for Ising-type systems this contribution is usually negligible [14]. 
Since our model presents a short-ranged interactions, 3D Ising criticality is 
expected. 

2.3 Freezing Line 

When a solid phase is involved, techniques with insertion and removal of parti- 
cles are highly inefficient since the insertion of a particle in a high density phase 
is strongly restricted. Hence, to locate the freezing line we carried out simula- 
tions in the canonical ensemble at high density. A system partially crystallized 
at the desired temperature and p = 0.90 is then expanded to lower density, 
by small steps {dp = 0.01) followed by an equilibration period. The crystal- 
lization degree is measured after equilibration for every density, by means of 
the global order parameter [15,16]: 



where the angular brackets denote an ensemble average. QemiT^lj) is calculated 

for the neighbors j of particle i: a neighbor is one particle within a given 
distance Vq = 1.5 from i. rlj is the vector which joins two neighboring particles 
and Qerra is related to the spherical harmonics: 



Note that Qq is invariant with respect to the reference coordinate system 
chosen. This global order parameter is close to zero in the fluid phase and has 
a non zero value when any degree of crystallization appears in the sample [15]. 
The freezing fine is then located as the lowest density with a noticiable degree 
of crystallization. 

2.4 Simulation Details 

We firstly study the gas-liquid coexistence using Gibbs Ensemble Monte Carlo 
(GEMC) simulations [17]. In the GEMC technique, two independent boxes are 
hold at the same temperature, pressure and chemical potential by allowing 
exchanges of particles and volume between the boxes, but keeping the total 
volume and particle number constant throughout the simulation. We compare 
GEMC results using N = 2000 particles and N = 432 to study the effects of 
finite size on the coexistence lines and the critical points. 




(2) 



(3) 
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Finite Size Scaling, based on the analysis of BW, have been performed by 
means of Grand Canonical Monte Carlo (GCMC) simulations aided with his- 
togram reweighting techniques. Phase diagram with L* = 10 were obtained 
for Ka — 6 , Ka — 10 and Ka — 15. Grand Canonical Simulations comprised 
25 • 10^ steps for L* = 8 and L* = 9 and 75 • 10^ steps for L* = 10 and 
L* — 12, each step consisted oi 2 < N > attempts to insert or remove a 
random particle. 

To estimate the freezing line we use standard NVT Monte Carlo simulations 
with N = 1024 particles. The particles are enclosed in a cubic box with 
periodic boundary conditions. 

3 Results and Discussion 

As we have discussed in previous papers [5,7], this model undergoes a gas- 
liquid transition in the low T-low p region; where high correlations between 
oppositely charged colloids arise, as in the RPM [6]. The critical temperature 
presents a non monotonic behaviour with the interaction range, what can be 
rationalized considering charge correlations [7] . On average, each charged col- 
loid is sorrounded by a layer of oppositely charged particles, followed by a 
layer of similarly charged colloids and so on. Therefore, when the salt con- 
centration is increased, the attractive term is shortly screened (because the 
opposite charged particles are in contact); contrary, the repulsive contribu- 
tions are strongly screened since longer distances separate unequal colloids. 
The system, thus, gains energy upon increasing na, resulting in an increase of 
the critical temperature (opposite to the behaviour of monocomponent sys- 
tems). At high enough k, the repulsive interactions are completely screened, Tc 
decreases as in a monocomponent system, since only the first layer of particles 
(with opposite charge) interacts, leaving a maximum at Ka ~ 10. An island 
of phase separation (gas-liquid) is then predicted for this model. 

As presented in the introduction, the aim of this paper is to study whether 
the triple temperature is lower than the critical one around Ka values which 
are close to the maximum Tc. We will study thus the liquid-gas transition, and 
then the freezing line. For the latter, we need first to form stable crystals at 
different temperatures. 

3.1 Liquid-gas transition 

GEMC and GCMC results for the liquid-gas transition at Ka = 6 and Ka = 
10 are presented in Fig. 1. Note that, contrary to monocomponent systems, 
the critical point is at higher temperature for the system with the shorter 
interaction range. To check for finite size effects, we compare in this figure 
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Fig. 1. Gas-liquid transition for Ka = 6 (filled symbols) and Ka = 10 (open symbols): 
GEMC results with N = 432 (squares) and N = 2000 (triangles) and GCMC results 
(diamods). The critical point in the GEMG simulations is calculated by means of 
the law of rectilinear diameters. 

simulations with = 432 and = 2000 particles (GEMC), and with those 
from GCMC simulations aided with reweighting techniques with L* = 10. 
Although the different estimations of the critical temperautures coincide quite 
well, bigger differences for the critical densities are noticed, in agreement with 
previous comparisons for the RPM and other models [18,19]. The coexistence 
curves coincide for the three cases far away from the critical point, and differ 
slightly only close to the critical point, where the correlation length becomes 
similar to (or even larger than) the simulation box. 

The behaviour of the critical parameters with the system size is studied in Fig. 
2 by means of GCMC simulations and using the corresponding scaling laws 
for the temperature (upper panel) and density (lower panel). Since the inter- 
actions are short ranged, the scaling exponents are taken from the 3D Ising 
model. This assumption is further supported by the analysis of the marginal 
distribution (see below). Both the critical temperature and density show only 
slight dependences on the system size, and the values extrapolated to an in- 
finite system differ very little from the critical values for the finite systems 
presented. 

In GCMC criticalicity is recognized when the marginal distribution of the 
simualtion coincides with that of the 3D Ising model, that is, the universality 
class must be assumed a priori. The best way of working would be to compare 
the distribution with different universal distributions, but for other models 
the critical distributions are not known. Therefore, we can only carry out a 
study of compatibility between our results and the 3D Ising universality class, 
which, nevertheless, is expected, due to the short range of the interactions. 
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Fig. 2. Critical temperature (upper panel) and density (lower panel) as a function 
of the system size, using the corresponding scaling laws and the exponents from the 
3D-Ising universality class. Black circles: Ka = 6, red squares: na = 10 and green 
diamonds: na = 15. 



The marginal distribution function {Pl[x)) is mapped onto the 3D Ising dis- 
tribution in the inset of Fig. 3 (lower panel), and the critical conditions are 
determined when the match between both distributions is achieved. Note that 
the distributions obtained from simulations match very well the universal 
Ising one. Additionally, using the relation ~ L^^^'^ the ratio between the 
critical exponents (3 and u can be obtained. The values computed are simi- 
lar to those of the Ising model {(3/v = 0.518) for the three ranges studied: 
13 /u = 0.523(14), 0.509(23) and 0.503(20) for Ka = 6,10 and 15, respectively 
(upper panel in Fig. 3). 

To finish the study of the GCMC data, we follow the analysis which was 
performed by Camp and Patey [20]. They showed that for systems with 3D 
Ising universality class, the relation = BL^/" should be fulfilled (where B 
is a system-dependent parameter). In the lower panel of Fig. 3, the results 



7 



J I I I I I L 



2 2.1 2.2 2.3 2.4 2.5 

In (L) 



CO 10- 




Fig. 3. Critical behavior study for Ka = 6 (circles), na = 10 (squares) and na = 15 
(diamods). Higher panel: estimation of the critical ratio j3/u. Lower panel: Am{L) 
versus L^/"". Blue lines are lineal fits for the numerical results. Inset: order param- 
eter probability distribution Pl{x) for = 10 and box sizes: L = 8 (circles), L = 9 
(squares), L = 10 (diamonds) and L = 12 (triangles); solid line for 3D Ising model. 



for this analysis are plotted. Note that for the three na values studied the 
results can be fitted with straight lines crossing the origin, indicating again 
the compatibility between our results and the Ising criticalicity. Because we 
have always found compatibility between our GCMC results and the Ising 
universality class, we can state that the critical parameters and the coexistence 
curve from GCMC are expected to be more accurate than the GEMC ones. 
However, since we are interested in the stability of the liquid phase, if the 
triple temperature is much lower than Tc, we can indistincly use GCMC or 
GEMC results because far away from the critical point both coexistence curves 
coincide (Fig. 1). 




Fig. 4. Snapshot of the system at T* = 0.17 and p* = 0.90 for kg = 6, showing 
liquid- crystal coexistence. Note the ordering of the particles, every particle has eight 
neighbours with opposite charge. This systems yields Qq = 0.41 

3.2 Structure of the crystal 

In the RPM, three different crystal phases have been reported [10]: disordered 
FCC at high temperatures, where the ions are randomly located; ordered 
BCC or CsCl in the low T regions, with each ion in the center of a cubic box 
surrounded by eight opposite ions in the vertexes of the cubic box, which, upon 
compresion suffers a first order transition to an ordered FCC structure. Due 
to the analogies between the colloidal mixture and the RPM, similar phases 
are to be expected here. Since we only attempt to nail down the freezing line, 
we will not discuss here the possibility of the transition between CsCl and 
ordered-FCC structures. 

In Fig. 4 we present a snapshot of the system at T* = 0.17 and p* = 0.90 
for Kcr = 6, which has very clearly crystallized into a CsCl structure. This 
temperature is very similar to the critical one for this range, which corre- 
spond to "low temperatures" , in the comparison with the RPM results. This 
picture compares nicely with experimental ones for symmetrical mixtures of 
charged colloids presented recently in similar conditions [8,9]. Similar crystal- 
lites were observed at other low temperatures and not-too-high density. At 
lower temperatures, however, the crystallite is in coexistence with a vapour 
phase, instead of a liquid phase, signalling the triple temperature. The system 
in the snapshot above and a system with T* = 0.13 and p* = 0.3 (gas-crystal 
coexistence) are analysed in Fig. 5. 

In the upper panel, the partial radial distribution functions are shown (black 

line for (r) and green line for (yf++(r)). Note that the system is composed of 

alternating layers of opposite charges, similarly to the structure of the liquid 
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Fig. 5. Upper panel: partial radial distribution functions for the system at T* = 0.13 
and p* = 0.30 for na = 6 (continuous lines) and at T* = 0.17 and p* = 0.90 

(broken lines), (r) with black line and g++ with green line. Lower panel: for both 

states, histograms of the number of particles inside a sphere centered in each particle 
with radius Rc: small cirlces for Rc = 1.15 (first layer of particles); squares for 
Rc = 1-15 considering only neighbors of oppositely charged; diamonds for Rc = 1.4 
(second layer of particles) and triangles for Rc = 1-85 (third layer of particles). 

phase. Noteworthly, the state with the lowest temperature yields narrower 
peaks, since the vapor is very dilute and it does not contribute noticeably to 
g{r), whereas the liquid smears out the crystal peaks at higher temperatures. 

The distribution of the number of particles inside a sphere centered in each 
particle with radius Rc is presented for both states in the lower panel of 
Fig. 5 (the central particle is not counted). When Rc = 1.15, only nearest 
neighbours are considered, i.e. the number of particles inside the first peak 
in g(r). Since the first and second layers overlap (see the upper panel), some 
colloids are surrounded by 10 particles. Nevertheless, when only particles with 
opposite sign are exclusively taken into account, the maximum number of 
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Fig. 6. Higher panel: total radial distribution function at T* = 2 and p* = 0.95 for 
Ka = 6 (fluid- crystal coexistence). Lower panel: as in Fig. 5: squares for Rq = 1.25 
and diamonds for Rq = 1.6. 

nearest neighbours is 8. That is, each sphere is surrounded at most by eight 
oppositely charged particles. Taking Rc = 1.4, the particles inside the first 
and second layers are being counted. Now the maximum number of particles 
inside such a sphere is 14. Finally, setting Rc = 1.85, we consider up to 
the third layer of particles and 24 particles were found. All these data are 
compatible with a BCC strucutre (as in the CsCl), but not with the FCC 
structure. As in the RPM, in the low T regions the fiuid (liquid or vapour) 
coexists with an ordered-BCC crystal. The stability of such crystals is due to 
the energetic terms which favor the formation of BCC-like crystals because the 
number of contacts with particles of opposite sign is higher than in the FCC 
structure, even though the entropic contribution aims for the more compact 
FCC structure. 

At higher temperatures, the results obtained are completely different. Fig. 
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6 shows the radial distribution function and the distribution of the number 
of particles inside spheres with different radii (T* = 2 and p* = 0.95). The 
partial structure is now random, that is, structuration in layers was not found. 
The distribution of particles is (lower panel in Fig. 6): 12 particles in the first 
layer and 6 particles in the second one, corresponding to the disordered-FCC 
crystal. Both structures shown here fully agree with the recent experimental 
results on the same system [8,9] and with the expectations from the RPM [10]. 

3.3 Stability of the liquid 

The freezing line can be determined using the procedure described in Section 
II: we take a spontaneuosly formed crystal seed at high enough density and 
decrease the density, with small steps, down to the point where the crystal 
is completely melted. Throughout this path, the global orientation parameter 
Qe was used as the order parameter. Fig. 7 shows the numerical estimations 
of the freezing line together with the gas-liquid coexistence points for three 
different kct values (ku = 6, where the critical temperature grows with k; 
na = 10, Tc reaches its maximum value and kct = 20, T^. behaves as in mono- 
component systems). For these three values studied here, the liquid phase is 
stable with respect crystallization as we can observe in Fig. 7 in a narrow 
window of temper aures. Interestingly, the freezing line behaves monotonically 
with Ka, moving to lower density (or higher temperature), contrary to the 
non-monotonous behaviour of the liquid-gas transition. 

The triple point is set by the crossing between the freezing line and the liquid 
branch. Fig. 8 plots the critical and triple temperatures for the colloidal sys- 
tem and its atomic analog, the ionic fluid. Note that the stable gas-liquid gap 
decreases as K.a increases. Extrapolating the present results, the triple temper- 
ature is expected to be equal to the critical one around Ka ~ 25, thus implying 
that the liquid-gas transition is mctastable with respect to crystallization for 
shorter interaction ranges. This prediction cannot be proved straight away 
because more sophisticated algorithms should be used to sample correctly the 
whole space phase for systems with such short-ranged interactions [21]. 



4 Conclusions 

We have studied the stabihty of the hquid phases for the colloidal analog of the 
ionic fluid, where the system was modeled by the effective DLVO interaction 
potential. This system presents richer behavior than the ionic one [5,7,8,9], 
because the range of the interaction can be tuned and the electroneutrality is 
imposed by construction. The liquid phases in these mixtures are stable for 
the range of the interaction studied, but the existence of the liquid is restricted 
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Fig. 7. Gas-liquid coexistence curves from GEMC ( open circles and filled circles for 
the critical points) and freezing line (filled triangles) for Ka = 6 (higher panel), 
K,a = 10 (middle panel) and Ka = 15 (lower panel). Note that the T*-scales differ 
for different panels. 



to narrower T-regions as the interaction range is decreased. 

BCC crystals have been found at low T, but at higher temperatures, FCC 
crystals appear, when the correlations between oppositely colloidal particles 
are not important. These results confirm that the phase diagram of a symmet- 
rical colloidal mixture presents similar aspect to that of the RPM (while the 
liquid is stable). Even, experimental studies support such a statement [8,9]. 
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Fig. 8. Squares: critical temperatures of the colloidal mixture from GEMC (triangles 
from GCMC) and the ionic fluid [19]. Stars: triple temperatures (estimated value 
for Ka = 6). The symbols at na = mark the RPM values [10]. 
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